Thierry Phénix & Ladislas Nalborczyk
21 mai 2019
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 : Modèle linéaire à 1 prédicteur continu
Cours 4 : Modèle linéaire multivarié
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Prédicteurs catégoriels et Interactions
Cours 9: Modèles multi-niveaux
Cours 10: Data Hackaton
- The first idea is that Bayesian inference is reallocation of credibility across possibilities.
- The second foundational idea is that the possibilities, over which we allocate credibility, are parameter values in meaningful mathematical models.
Kruschke (2015)
3 configurations possibles:
\( \longrightarrow~ \) Présenter le principe de base de l'échantillonnage : les MCMC
\( \longrightarrow~ \) Présenter les méthodes classiques (point de vue historique)
\( \longrightarrow~ \) Montrer les forces mais aussi les faiblesses de ces méthodes
\( \longrightarrow~ \) Donner des outils de contrôle sur ces méthodes
\( \longrightarrow~ \) Appliquer ces méthodes à un cas simple
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 : Modèle linéaire à 1 prédicteur continu
Cours 4 : Modèle linéaire multivarié
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Prédicteurs catégoriels et Interactions
Cours 9: Modèles multi-niveaux
Cours 10: Data Hackaton
Markov chain Monte Carlo
\( \longrightarrow~ \) Échantillonnage aléatoire
\( \longrightarrow~ \) Le résultat est un ensemble de valeurs du paramètres
Markov chain Monte Carlo
\( \longrightarrow~ \) Les valeurs sont générées sous forme de séquence (liaison de dépendance)
\( \longrightarrow~ \) Indice temporel pour identifier la place dans la chaine
\( \longrightarrow~ \) Le résultat est de la forme \( ~:~\theta^1, \theta^2,\theta^3, \cdots, \theta^t \)
Markov chain Monte Carlo
\( \longrightarrow~ \) La valeur de paramètre générée ne dépend que de la valeur du paramètre précédent
\( \longrightarrow~P(\theta^{t+1}~|~\theta^{t},\theta^{t-1}, \cdots,\theta^{1})~=~P(\theta^{t+1}~|~\theta^{t}) \)
Le terme méthode de Monte-Carlo, désigne une famille d'algorithmiques visant à calculer une valeur numérique approchée en utilisant des procédés aléatoires, c'est-à-dire des techniques probabilistes.
Wikipedia
Cette méthode a été inventé en 1947 par Nicholas Metropolis, et publié pour la première fois en 1949 dans un article coécrit avec Stanislaw Ulam.
Les méthodes Monte Carlo designe une famille d'algorithmes qui ont pour but de calculer des valeurs numériques approchées à partir de procédés aléatoires
\( \longrightarrow~ \) Utiliser une méthode Monte Carlo pour calculer une valeur approchée de \( \pi \)
La probabilité que le point \( M \) appartienne au disque est \( \pi⁄4 \).
Les méthodes Monte Carlo designe une famille d'algorithmes qui ont pour but de calculer des valeurs numériques approchées à partir de procédés aléatoires
\( \longrightarrow~ \) Utiliser une méthode Monte Carlo pour chercher un maximum (Optimisation)
Par exemple déterminer le maximum d'une fonction.
Les méthodes Monte Carlo designe une famille d'algorithmes qui ont pour but de calculer des valeurs numériques approchées à partir de procédés aléatoires
\( \longrightarrow~ \) Utiliser une méthode Monte Carlo pour approximer une distribution de probabilité
On connait les priors \( ~P(\theta_1) \) et \( ~P(\theta_2) \)
On connait la likelihood \( ~P(D~|~\theta_1, \theta_2) \)
Mais on ne sait pas calculer la distribution postéreure \( ~P(\theta_1, \theta_2~|~D)=\frac{P(D~|~\theta_1, \theta_2)~P(\theta_1)P(\theta_1)~}{P(D)} \)
DONC on sait calculer la distribution postérieure à une constante près car \( P(D) \) est une constante.
On va évaluer la distribution postérieure à partir de notre connaissance sur les priors et la likelihood avec une méthode d'échantillonnage
Considérons un exemple simple
Soit un paramètre \( \theta \) qui ne prend que 7 valeurs, et une répartition de ce paramètre :
Question : Quelle est la distribution de \( \theta~ \)?
Remarque : mais c'est évident… ICI
Estimation de cette distribution par tirage aléatoire
On tire aléatoirement un grand nombre de point au hasard parmi ces 28 cases!
(Comme pour le calcul de \( \pi \))
Le résultat du tirage \( (4,6,2,3,7,7,5,3,6,5,1, \dots) \)
Que peut-on dire de cette méthode dans le cadre d'une distribution :
\( ~\bullet~ \) Elle converge vers la distribution
\( ~\bullet~ \) Aucun contrôle sur la vitesse de convergence
\( ~\bullet~ \) Nécessite beaucoup de points (7 valeurs de paramètres et 100 tirages)
\( ~\rightarrow~ \) En l'état, résultat pas concluant !
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 : Modèle linéaire à 1 prédicteur continu
Cours 4 : Modèle linéaire multivarié
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Prédicteurs catégoriels et Interactions
Cours 9: Modèles multi-niveaux
Cours 10: Data Hackaton
Cet algorithme a été présenté pour la première fois en 1953
Les auteurs :
Nicholas Metropolis, Arianna W. Rosenbluth, Marshall Rosenbluth, Augusta H. Teller, et Edward Teller
Marshall Nicholas Rosenbluth
(5 February 1927 – 28 September 2003)
Physicien Americain (spécialiste du plasma)
Première méthode
Le problème n'est pas la convergence, mais la vitesse à laquelle la méthode converge
Pour augmenter la vitesse de convergence, il faudrait faciliter l'accès aux valeurs de paramètres les plus représentées.
Principe
\( \bullet~ \) On fait une proposition de déplacement sur la base de la valeur courante du paramètre
\( \bullet~ \) on réalise un tirage aléatoire pour accepter ou rejeter la nouvelle position
2 Idées
\( \bullet~ \) la proposition doit favoriser les valeurs de paramètre les plus probables
\( ~~\longrightarrow~ \) On parcourt plus souvent ces valeurs de paramètres
\( \bullet~ \) la proposition doit se limiter aux valeurs adjacentes du paramètre courant
\( ~~\longrightarrow~ \) On augmente la vitesse de convergence en restant là où se trouve l'information
\( \longrightarrow~ \) Sélectionner un point de départ
On peut sélectionner n'importe quelle valeur
\( \longrightarrow~ \) Sélectionner un point de départ
On peut sélectionner n'importe quelle valeur
\( \longrightarrow~ \) Faire une proposition de déplacement
Faire un tirage et proposer un déplacement
\( \longrightarrow~ \) Sélectionner un point de départ
On peut sélectionner n'importe quelle valeur
\( \longrightarrow~ \) Faire une proposition de déplacement
Faire un tirage et proposer un déplacement
\( \longrightarrow~ \) Accepter ou rejeter la proposition de déplacement
Faire un tirage suivant la probabilité :
\[ P_{déplacement}=min \left(\frac{P(\theta_{proposé})}{P(\theta_{courant})},1 \right) \]
Remarque : Plus la valeur de \( \theta \) est grande et plus petite est la probabilité de rester sur place
\( \longrightarrow~ \) Sélectionner un point de départ
On peut sélectionner n'importe quelle valeur
\( \longrightarrow~ \) Faire une proposition de déplacement
Faire un tirage et proposer un déplacement
\( \longrightarrow~ \) Accepter ou rejeter la proposition de déplacement
Faire un tirage suivant la quasi probabilité :
\[ P_{déplacement}=min \left(\frac{P(\theta_{proposé})}{P(\theta_{courant})},1 \right) \]
\( \longrightarrow~ \) La position calculée devient la nouvelle position de départ
\( ~\bullet~ \) La likelihood est donnée par : \( \color{orangered}{P(z~|~\theta,N)=\theta^z~(1-\theta)^{(N-z)}} \)
\( ~\bullet~ \) Le prior est donné par : \( \color{steelblue}{P(\theta~|~a,b) \propto \theta^{(a-1)}(1-\theta)^{(b-1)}} \)
\( ~\bullet~ \) Le paramètre que l'on cherche à estimer \( \theta \) prend ses valeurs dans l'intervalle \( \left[0,1\right] \)
Problème 1: Comment définir la proposition de déplacement?
Le déplacement est modélisé par une distribution normale :\( ~\Delta\theta~\sim ~\mathcal{N}0,\sigma) \)
\( \longrightarrow~ \) la moyenne \( \mu \) vaut \( 0 \) : le déplacement se fait autour de la valeur courante du paramètre
\( \longrightarrow~ \) La variance reste à déterminer, elle contrôle l'éloignement de la nouvelle valeur
Problème 2: Quelle probabilité utiliser pour accepter ou refuser le déplacement?
Nous utilisons le produit de la likelihood et du prior : \( P(\theta)~\propto~ \color{orangered}{\theta^z~(1-\theta)^{(N-z)}}\color{steelblue}{\theta^{(a-1)}(1-\theta)^{(b-1)}} \)
La probabilité d'accepter le déplacement : \( ~P_{déplacement}=min \left(\frac{P(\theta_{cur}+\Delta\theta)}{P(\theta_{cur})},1 \right) \)
REMARQUE : le rapport \( \frac{P(\theta_{cur}+\Delta\theta)}{P(\theta_{cur})} \) est le même que l'on utilise la distribution postérieure ou le produit prior par likelihood !
\( \longrightarrow~ \) Sélectionner un point de départ
\( ~\bullet~ \) Il faut choisir \( ~\theta \in \left[0,1\right] \)
\( ~\bullet~ \) Seule contrainte \( ~P(\theta_{init}) \ne 0 \)
\( \longrightarrow~ \) Faire une proposition de déplacement
\( ~\bullet~ \) Faire un tirage suivant \( ~\mathcal{N}0,\sigma) \)
\( \longrightarrow~ \) Accepter ou rejeter la proposition de déplacement
Faire un tirage suivant la probabilité :
\[ P_{déplacement}=min \left(\frac{P(\theta_{cur}+\Delta\theta)}{P(\theta_{cur})},1 \right) \]
\( \longrightarrow~ \) La position calculée devient la nouvelle position de départ
Deux indice permettre d'évaluer la qualité de l'échantillonnage :
\( \rightarrow~ \) Le rapport entre le nombre de déplacements proposés et le nombre de déplacements accéptés
\( \rightarrow~ \) L'effective sample size
C'est le nombre de déplacement qui ne sont pas corrélés avec les précédents
\( \rightarrow~ \) Toutes les propositions de déplacement (ou presque) sont acceptées
\( \rightarrow~ \) Peu de valeurs effectives
Il va falloir beaucoup d'itérations pour avoir un résultat satisfaisant
\( \rightarrow~ \) Les propositions de déplacement sont rarement acceptées
\( \rightarrow~ \) Peu de valeurs effectives
Il va falloir beaucoup d'itérations pour avoir un résultat satisfaisant
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 : Modèle linéaire à 1 prédicteur continu
Cours 4 : Modèle linéaire multivarié
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Prédicteurs catégoriels et Interactions
Cours 9: Modèles multi-niveaux
Cours 10: Data Hackaton
Cet algorithme a été présenté pour la première fois en 1984
par les frêres Geman
L'algorithme Gibbs pour évaluer une distribution posterieure ne fonctionne qu'avec des modèles composés d'au moins deux variables.
On considère une situation dans laquelle on dispose de deux pièces de monaie.
A chaque pièce on associe un paramètre ou biais \( ~\theta_1 \) et \( ~\theta_2 \) (la probabilité d'obtenir \( Face \)).
On cherche à estimer les biais de ces deux pièces : \( ~P(\theta_1, \theta_2~|~D ) \)
Nos variables \( \theta_1 \) et \( \theta_2 \) sont indépendantes et identiquement distribuées (iid)
On considère que chacun des biais est décrit par une distribution beta :
\( \theta_1 \sim beta(\theta_1~|~a_1,b_1) \) et \( ~\theta_2 \sim beta(\theta_2~|~a_2,b_2) \)
On peut donc calculer la distribution a priori de \( ~P(\theta_1, \theta_2) \) :
\[ \begin{align} P(\theta_1, \theta_2) &= P(\theta_1)~P(\theta_2) \\ &= \theta_1^{(a_1-1)}~(1-\theta_1)^{(b_1-1)}\theta_2^{(a_2-1)}~(1-\theta_2)^{(b_2-1)} \end{align} \]
On dispose d'un ensemble de lancers pour les deux pièces
On appelle \( ~D_1 \) les données issues de la première pièce, composé de \( ~z_1 \) Face sur \( ~N_1 \) lancers.
On appelle \( ~D_2 \) les données issues de la seconde pièce, composé de \( ~z_2 \) Face sur \( ~N_2 \) lancers.
On note les données : \( D = \{ z_1, N_1, z_2, N_2 \} \)
\[ \begin{align} P(D~|~\theta_1,\theta_2) &= \prod_{y_{1i} \in D_1} P(y_{1i}~|~\color{orangered}{\theta_1})~ \prod_{y_{2i} \in D_2} P(y_{2i}~|~\color{orangered}{\theta_2}) &&\mbox{Règle d'indépendance}\\ &= \theta_1^{z_1}~(1-\theta_1)^{(N1-z_1)}\theta_2^{z_2}~(1-\theta_2)^{(N2-z_2)}~ \end{align} \]
\[
\begin{align}
\color{purple}{P(\theta_1,\theta_2~|~D)} &= \color{orangered}{P(D~|~\theta_1,\theta_2)}~\color{steelBlue}{P(\theta_1,\theta_2)}/\color{green}{P(D)} \\
&= \color{orangered}{\theta_1^{z_1}~(1-\theta_1)^{(N1-z_1)}\theta_2^{z_2}~(1-\theta_2)^{(N2-z_2)}}~\color{steelBlue}{P(\theta_1,\theta_2)}/\color{green}{P(D)} \\
&= \frac{\theta_1^{(\color{orangered}{z_1}+\color{steelBlue}{a_1-1})}(1-\theta_1)^{(\color{orangered}{N1-z_1}+\color{steelBlue}{b_1-1})}\theta_2^{(\color{orangered}{z_2}+\color{steelBlue}{a_2-1})}(1-\theta_2)^{(\color{orangered}{N2-z_2}+\color{steelBlue}{b_2-1})}~}{\color{green}{P(D)}~\color{steelBlue}{B(a_1,b_1)~B(a_2,b_2)}}
\end{align}
\]
Donc finalement
\[
\begin{align}
\color{purple}{P(\theta_1,\theta_2~|~D)} =& \color{purple}{beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1)} \\
& \color{purple}{beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)}
\end{align}
\]
On voit que lorsque le prior est le produit de deux distributions beta indépendantes, la distribution postérieure est le produit de deux distributions beta indépendantes
\( ~\longrightarrow~ \) La position suivante peut-être n'importe où sur la grille
On voit que l'efficacité de l'algorithme Metropolis repose sur la proposition de déplacement.
Cette fonction peut-être difficile à calibrer:
\( \longrightarrow~ \) on voudrait que sd soit petit dans les zones à forte densité
\( \longrightarrow~ \) au contraire, on voudrait que sd soit grand dans les zones à faible densité
L'algorithme Gibbs résout ce problème
L'algorithme de Gibbs parcourt l'espace des paramètres comme l'algorithme Metropolis
Mais plutôt que de modifier tous les paramètres en même temps, il ne modifie qu'un paramètre après l'autre
\( \longrightarrow~ \) la proposition de déplacement est basé sur la fonction de densité conditionnelle
\( \longrightarrow~ \) De ce fait, le déplacement est toujours accepté
\( \longrightarrow~ \) Sélectionner un point de départ \( (\theta_{1init},\theta_{2init}) \)
On peut sélectionner n'importe quelle valeur
\( \longrightarrow~ \) Sélectionner le premier paramètre par exemple \( \theta_1 \)
\( \longrightarrow~ \) Générer une valeur aléatoire de \( \theta_1 \)
Faire un tirage directement suivant la probabilité \( P(\theta_1~|~\theta2,D) \)
\( \longrightarrow~ \) Générer une valeur aléatoire de \( \theta_2 \)
Faire un tirage suivant la probabilité \( P(\theta_2~|~\theta1,D) \)
Avec pour valeur de \( \theta_1 \) la valeur calculée précédemment
\( \longrightarrow~ \) On boucle sur les deux dernières étapes pour parcourir
Pour utiliser l'algorithme, il faut pouvoir calculer toutes les fonctions de densité conditionnelle
Dans notre exemple :
\[ \begin{align} P(\theta_1~|~\theta_2, D) &= \frac{P(\theta_1,\theta_2~|~D)}{\int P(\theta_1,\theta_2~|~D)~\mathrm d\theta_1} \\ &= \frac{beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1)~beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)}{\int beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1)~beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)~\mathrm d\theta_1} \\ &= \frac{beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1)~beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)}{beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)~\color{orangered}{\int beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1)~\mathrm d\theta_1}} \\ &= \frac{beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)}{beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)}~beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1) \\ &= beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1) \end{align} \]
En conclusion,
\[ P(\theta_1~|~\theta_2, D) = P(\theta_1~|~D) \]
L'algorithme Gibbs est particulièrement utile lorsque la distribution jointe \( P(\{\theta_i\}~|~D) \) n'est pas calculable ou facile à échantillonner, mais que les densité conditionnelle \( P(\theta_i~|~\{\theta_{j \ne i}\}, D) \) le sont.
L'algorithme Gibbs peut-être vu comme un cas spécial de l'algorithme Metropolis pour lequel la proposition de saut dépend de la position dans l'espace des paramètres et du paramètre sélectionné.
En général, L'algorithme Gibbs améliore les performances de l'algorithme Metropolis
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 : Modèle linéaire à 1 prédicteur continu
Cours 4 : Modèle linéaire multivarié
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Prédicteurs catégoriels et Interactions
Cours 9: Modèles multi-niveaux
Cours 10: Data Hackaton
Cet algorithme a été présenté pour la première fois en 1987
les auteurs :
Simon Duane, A.D. Kennedy, Brian Pendleton and Duncan Roweth
On a vu que l'algoritme Gibbs améliore les performances de convergence de l'algorithme Metropolis.
Mais il nécessite de pouvoir calculer les densités conditionnelles \( P(\theta_i~|~\theta_{j \ne i}, D) \)
L'algoritme est performant lorsqu'il y a indépendance des paramètres \( P(\theta_i~|~\theta_{j \ne i}, D) = P(\theta_i~|~D) \) Ce n'est pas le cas en général!!!
Pour le HMC, abandon de cette stratégie
L'algorithme Hamiltonian Monte Carlo est une amélioration de Metroplolis.
Une version particulière de l'algorithme Metropolis :
\( \longrightarrow~ \) la proposition de déplacement est toujours une gaussienne multivariée centrée sur la position courante
\( \longrightarrow~ \) la proposition de déplacement est une gaussienne qui garde toujours la même forme (variance constante)
L'adaptation de la proposition de déplacement en fonction de la position courante du paramètre.
\( \longrightarrow~ \) Sélectionner un point de départ
\( ~\bullet~ \) Il faut choisir \( ~\theta \in \left[-3,3\right] \) dans cet exemple
\( \longrightarrow~ \) Il faut définir une fonction de potentiel
On prend l'opposé du logarithme de la fonction de densité
\( \longrightarrow~ \) On peut générer un déplacement
On présente ici un ensemble de point pour voir la tendance
\( \longrightarrow~ \) L'histograme présente la tendance
On voit que le mode de la distribution de la future position n'est pas aligné avec la position du point de départ
\( \longrightarrow~ \) Sélectionner un point de départ
On peut sélectionner n'importe quelle valeur de \( \theta \) dans l'espace des paramètres
\( \longrightarrow~ \) Générer une nouvelle position à partir de la fonction de potentiel
\( ~\bullet~ \) L'algorithme génère une proposition par analogie au lancé d'une bille sur la fonction de potentiel
\( ~\bullet~ \) Le pas de déplacement, ainsi que le nombre de déplacement est fixé à l'avance
\( ~\bullet~ \) La force avec laquelle on lance la bille est déterminée par une loi normale centrée
\( \longrightarrow~ \) Accepter ou rejeter la proposition de déplacement
Faire un tirage suivant la probabilité suivante en considérant \( \phi \) comme le moment associé à la bille
\[ P_{déplacement}=min \left(\frac{P(\theta_{proposé}~|~D)~P(\phi_{proposé})}{P(\theta_{courant}~|~D)~P(\phi_{courant})},1 \right) \]
\( \longrightarrow~ \) on enregistre la nouvelle position et on recommence…
La faiblesse est liée à la génération de la proposition de déplacement.
quand le pas de déplacement est trop petit
ou
quand le nombre de pas est trop grand
La faiblesse est liée à la génération de la proposition de déplacement.
L'influence du premier potentiel choisi
- L'algorithme HMC est particulièrement utile lorsque le modèle contient un grand nombre de paramètres et quand les paramètres sont corrélés entre eux.
- Il nécessite un paramétrage assez fin. Des logiciels proposent des méthodes clés en main très efficaces.
Exemple : STAN
- Le package brms (déjà vu aux chapitres précédents) utilise cet algorithme … \( \longrightarrow~~~ \) fonction : brm()
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 : Modèle linéaire à 1 prédicteur continu
Cours 4 : Modèle linéaire multivarié
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Prédicteurs catégoriels et Interactions
Cours 9: Modèles multi-niveaux
Cours 10: Data Hackaton
A. Le temps de calcul limité !
B. ces méthodes nécessitent un paramétrage fin
\( \longrightarrow~ \) Algorithme Gibbs
\( ~\bullet~ \) la variance de la distribution normale de la proposition
\( \longrightarrow~ \) Algoorithme HMC
C'est le paramétrage du calcul de la position suivante qui est délicat :
\( ~\bullet~ \) le pas
\( ~\bullet~ \) la longueur du pas
\( ~\bullet~ \) la force appliquée à la bille
Ces méthodes produisent des chaine de valeurs de paramètres (échantillons).
Chaque chaine produite doit respecter trois critères pour être validée
Elle doit-être représentative de la distribution postérieure
\( ~\bullet~ \) elle ne doivent pas dépendre du point de départ
\( ~\bullet~ \) elle ne doit pas être cantonnée à une région de l'espace des paramètres
Elle doit-être suffisamment longue pour assurer la précision et la stabilité du résultat
\( ~\bullet~ \) La tendance centrale et le HDPI calculés à partir de la chaine ne doivent pas changer si on relance la procédure
Elle doit-être générée de manière efficace
\( \longrightarrow~ \) Vérification visuelle des trajectoires
\( ~\bullet~ \) Les chaines doivent occuper le même espace
\( ~\bullet~ \) la convergence ne dépend pas du point de départ
\( ~\bullet~ \) Aucune chaine ne doit avoir de trajectoire particulière (cyclique)
\( \longrightarrow~ \) Vérification visuelle des densités
\( ~\bullet~ \) Les densités doivent se superposer
\( ~\bullet~ \) Les densités doivent-être régulières
Cette affichage ne montre que les 500 premières itérations
Les trajectoires ne se superposent pas au début (zone orange)
La densité est également affectée
C'est la Burn-in period
En pratique on supprime ces premières itérations
\( \longrightarrow~ \) Vérification numérique des chaines
le Shrink factor (en bas à gauche)
Correspond au “Rhat” de brms
C'est le rapport entre la variance inter chaines et intra chaines.
\( \bullet~ \) intuitivement la valeur doit tendre vers 1.0
\( \bullet~ \) en pratique elle est acceptable jusqu'à 1.1
Plus la chaine est longue et plus le résultat sera précis et stable
Si la chaine s'attarde sur chaque position, et que le nombre d'itrations reste le même, alors on perd en précision
Il lui faudra plus d'itérations pour arriver au même niveau de précision
\( \longrightarrow~ \) Mesure l'autocorrélation
L'autocorrelation est la corrélation de la chaine avec elle-même mais décalé de \( k \) itérations (lag)
\( \longrightarrow~ \) L'autocorrelation doit tendre vers 0 lorsque le lag augmente
La fonction d'autocorrélation est représentée pour chaque chaine (en haut à droite)
Un autre résultat rend compte de la précision de l'échantillon : l'effective sample size ESS
\[ ESS~=~\frac{N}{1+2\sum_k ACF(k)} \]
ESS représente la taille d'un échantillon non-autocorrélé extrait de la somme de toutes les chaines
Pour une précision raisonnable du HDI, il est recommandé d'avoir un ESS autour de 10000
L'erreur standard d'un échantillon est donné par le calul : \( SE~=~SD/ \sqrt{N} \)
plus N augmente plus l'erreur standard diminue
On peut généraliser cette idée aux chaines de markov
\[ MCSE~=~SD/ \sqrt{ESS} \]
Pour une précision raisonnable de la tendance centrale, il faut que cette valeur soit petite
Dans l'exemple, l'ESS est petite mais malgré ça le MCSE est petite et le mode est assez stable
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 : Modèle linéaire à 1 prédicteur continu
Cours 4 : Modèle linéaire multivarié
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Prédicteurs catégoriels et Interactions
Cours 9: Modèles multi-niveaux
Cours 10: Data Hackaton
set.seed(8)
y <- rnorm(100, mean = 0, sd = 1)
b5.1 <-
brm(data = list(y = y,
intercept_1 = 1,
intercept_2 = 1),
family = gaussian,
y ~ 0 + intercept_1 + intercept_2,
prior = c(prior(uniform(-1e10, 1e10), class = b),
prior(cauchy(0, 1), class = sigma)),
inits = list(list(intercept_1 = 0, intercept_2 = 0, sigma = 1),
list(intercept_1 = 0, intercept_2 = 0, sigma = 1)),
iter = 4000, warmup = 1000, chains = 2,
seed = 8)
set.seed(8)
y <- rnorm(100, mean = 0, sd = 1)
b5.2 <-
brm(data = list(y = y,
intercept_1 = 1,
intercept_2 = 1),
family = gaussian,
y ~ 0 + intercept_1 + intercept_2,
prior = c(prior(normal(0, 10), class = b),
prior(cauchy(0, 1), class = sigma)),
inits = list(list(intercept_1 = 0, intercept_2 = 0, sigma = 1),
list(intercept_1 = 0, intercept_2 = 0, sigma = 1)),
iter = 4000, warmup = 1000, chains = 2,
seed = 8)
Remplacement du prior uniform par un prior normal